#This script relies on the Zelig package.  In particular, I use version 3.5.5 to produce these analyses.
#To install this version of Zelig: 
#Download 3.5.5 from here: http://cran.r-project.org/src/contrib/Archive/Zelig/ (into R's working directory).
# Then type this in R: 
#install.packages("http://cran.r-project.org/src/contrib/Archive/Zelig/Zelig_ 3.5.5.tar.gz", repos=NULL, type="source")
#Note you may need to delete your current installation of Zelig manually in the folder where it is located

rm(list=ls(all=TRUE))
library(readstata13)
library(Zelig)
library(car)
library(apsrtable)
library(readxl)
library(plyr)
library(xtable)
options(digits=2)
options(scipen = 10)

#SET WORKING DIRECTORY HERE
#setwd('~/Dropbox/Panel Effects/Race cases/Death Penalty and Race/Paper files/JSJ_submission/Replication_files')

###############################################################
#Data on number of death sentences over time
###############################################################
sentence.data <- read_excel('death_sentences_by_year.xlsx')

#MAKE FIGURE 1

pdf("death_sentences_over_time.pdf", width=10, height=6)

par(mfrow=c(1,1), mar=c(2,4.5,1,1))
plot(sentence.data$year, sentence.data$n.sentences, type="n", axes=F, xlab="", ylab="", lwd=1.5)
polygon (x = c(1973, 1973, 1995, 1995), y=par()$usr[c(3,4,4,3)],  col=gray(.85))
points(sentence.data$year, sentence.data$n.sentences, type="l")
axis(1, at=seq(1975,2020,5), cex.axis=1.4, mgp=c(2,.5,0))
axis(2, cex.axis=1.4, mgp=c(2,.5,0), las=1)
mtext("Number of death sentences per year", 2, srt=90, cex=1.5 ,line=3)
box()

dev.off()

###############################################################
#Public opinion data on death penalty
###############################################################

#GSS
gss <- read_excel('gss_death_penalty_by_race_for_analysis.xlsx')
head(gss)

pdf("gss_public_opinion.pdf", width=10, height=6)

par(mfrow=c(1,1), mar=c(2,4.5,1,1))
plot(gss$year, gss$overall, type="n", axes=F, ylim=c(30,80), xlab="", ylab="")
#shade 1973 to 1995
foo <- c( 1973, 1973, 1995, 1995) 

polygon (x = foo, y=par()$usr[c(3,4,4,3)],  col=gray(.85))
points(gss$year, gss$overall,type="l", col="dark green", lwd=1.5)
points(gss$year, gss$white,type="l", lty=2, col="red", lwd=1.5)
points(gss$year, gss$black,type="l", lty=6, col="black", lwd=1.5)
axis(1, at =seq(1975,2015,5), mgp=c(2,.5,0), cex.axis=1.4)
axis(2, at = seq(30,80,10),   cex.axis=1.4, mgp=c(2,.5,0), las=1)
mtext("Percent who favor death penalty", 2, srt=90, cex=1.5 ,line=3)
text(2005, 73.4, "Whites", cex=2, col="red")
text(2005, 63, "Overall", cex=2, col="dark green")
text(2010, 47, "Blacks", cex=2, col="black")

box()

dev.off()


############################################################################################################
#LIEBMAN CASE DATA FROM HERE ON OUT 
############################################################################################################
lieb.cd <- read.dta13("liebman_case_data_for_replication_script.dta", convert.underscore = TRUE)

#DESCRIPTIVE STATS In text

#how many cases in data?
nrow(lieb.cd)

#mean of cases where relief is granted
mean(lieb.cd$panel.relief)

#proportion of mixed panels
mean(lieb.cd$mixed.panel)

#hispanic panels and defendant
with(lieb.cd, mean(hispanic1 + hispanic2 + hispanic3>=1))
mean(lieb.cd$def.race=="hispanic")

##proportions of defendant race
prop.table(table(lieb.cd$def.race))
mean(lieb.cd$def.race=="black")
mean(lieb.cd$def.race=="white")

#victim distribution
table(lieb.cd$vic1.race)
table(lieb.cd$vic1.race)[2]/(length((lieb.cd$vic1.race)))#proportion of black victims

#where do mixed panels occur?
table(lieb.cd$mixed.panel, lieb.cd$state)
table(lieb.cd$mixed.panel, lieb.cd$circuit)

################
#TABLE 1: balance checks: check balance across cases with and without black judges
  #note: I omit panel liberalism and circuits from the published version of the table
#choose vars and remove observations with missing values
complete.lieb.cd <-  subset(lieb.cd, select=c("mixed.panel","black.def","state.evid.hearing","out.state.representation",
	"aggrav.mit.index", "offender.victim.index","fed.evid.hearing","number.claims",
	"year.sentence","panel.liberalism","number.dems", "circuit3", "circuit4", "circuit5", 
	"circuit6", "circuit7", "circuit8", "circuit9","circuit10", "circuit11")) 
complete.lieb.cd <- complete.lieb.cd[complete.cases(complete.lieb.cd),]

vars <-  names(complete.lieb.cd[names(complete.lieb.cd)!="mixed.panel"])
ttest.results <-lapply(complete.lieb.cd [,vars], function(x) t.test(x ~ complete.lieb.cd$mixed.panel, var.equal = TRUE))
ttest.results
#extract from list
var.names <- names(ttest.results)
var.names <- c("Black defendant", "State evidentiary hearing", "Out of state representation",
	"Aggravation/Mitigation index","Offender victim index","Federal evidentiary hearing",  
	"Number of claims", "Year sentenced", 
	"Panel liberalism",
	"Number of Democrats on panel", "3rd circuit", "4th circuit",
	  "5th circuit","6th circuit",  "7th circuit", "8th circuit", "9th circuit", 
		"10th circuit", "11th circuit")
mean.nonblack.panel <- sapply(ttest.results, "[[", "estimate" )[1,]
mean.mixed.panel <- sapply(ttest.results, "[[", "estimate" )[2,]
p.values <- round(sapply(ttest.results, "[[", "p.value" ),2)
#do multiple comparisons test on everything but number.dems
multiple.p <- round(p.adjust(p.values,c("BH")),2)
match.test <- data.frame(var.names, mean.nonblack.panel, mean.mixed.panel, p.values, multiple.p )
colnames(match.test) <- c("Variable", "Mean, non-mixed panels", "Mean, mixed panels", "p-value", "Corrected p-value")
rownames(match.test) <- NULL
tab <- xtable::xtable(match.test, label="balance_test", caption="Balance statistics for panels with and without black judges. 
	The	`p-value' column depicts the results of a t-test comparing the means relevant variable across cases heard by mixed and non-mixed panels. 
	The column `Corrected p-value' uses the method of citet{benjamini:1995} to adjust for multiple comparisons.
	After accounting for multiple comparisons, the only statistically significant difference is that panels with black judges feature more Democratic judges (and thus are more liberal overall), since
	all black judges in the data are Democrats.", size="footnotesize")

########################################################################################
#REGRESSIONS FOR TABLE 2
########################################################################################
#MAIN REGRESSIONS
#1) Just mixed.panel and black.def
	#Liebman
lieb.main.no.inter <- glm(panel.relief ~
		 panel.liberalism + state.evid.hearing + out.state.representation  + aggrav.mit.index  + offender.victim.index 
		+ fed.evid.hearing  + number.claims + year.sentence.trend
			+ circuit4 + circuit5 + circuit8 + circuit10 + circuit11
			 + black.def	+ mixed.panel
				,	data = lieb.cd,family=binomial(link="logit"))
summary(lieb.main.no.inter)		

#add  Interaction--raw model with no controls
	#LIEBMAN
lieb.main.with.inter.raw <- glm(panel.relief ~
		+ mixed.panel	+ black.def	+ mixed.panel:black.def
						,	data = lieb.cd,family=binomial(link="logit"))
summary(lieb.main.with.inter.raw)
lieb.main.with.inter.raw.b1b3 <- coef(lieb.main.with.inter.raw)["mixed.panel"] + coef(lieb.main.with.inter.raw)["mixed.panel:black.def"]
lieb.main.with.inter.raw.b1b3
lieb.main.with.inter.raw.hyp <- linearHypothesis(lieb.main.with.inter.raw, "mixed.panel + mixed.panel:black.def=0")
lieb.main.with.inter.raw.hyp

#now add interaction with black defendant, no liberalism
#LIEBMAN
lieb.main.with.inter <- glm(panel.relief ~
		+ mixed.panel	+ black.def	+ mixed.panel:black.def
		+ state.evid.hearing + out.state.representation  + aggrav.mit.index  + offender.victim.index 
		+ fed.evid.hearing  + number.claims +  year.sentence.trend
		+ circuit4 + circuit5 + circuit8 + circuit10 + circuit11 
				,	data = lieb.cd,family=binomial(link="logit"))
print(summary(lieb.main.with.inter, detail=T),digits=3)
lieb.main.with.inter.b1b3 <- coef(lieb.main.with.inter)["mixed.panel"] + coef(lieb.main.with.inter)["mixed.panel:black.def"]
lieb.main.with.inter.hyp <- linearHypothesis(lieb.main.with.inter, "mixed.panel + mixed.panel:black.def=0")
lieb.main.with.inter.hyp 

#now add liberalism
#LIEBMAN
lieb.main.with.inter.and.lib <- glm(panel.relief ~
		+ mixed.panel	+ black.def	+ mixed.panel:black.def + panel.liberalism
		+ state.evid.hearing + out.state.representation  + aggrav.mit.index  + offender.victim.index 
		+ fed.evid.hearing  + number.claims + year.sentence.trend
		+ circuit4 + circuit5 + circuit8 + circuit10 + circuit11
				,	data = lieb.cd,family=binomial(link="logit"))
print(summary(lieb.main.with.inter.and.lib, detail=T),digits=3)
lieb.main.with.inter.and.lib.b1b3 <- coef(lieb.main.with.inter.and.lib)["mixed.panel"] + coef(lieb.main.with.inter.and.lib)["mixed.panel:black.def"]
lieb.main.with.inter.and.lib.hyp <- linearHypothesis(lieb.main.with.inter.and.lib, "mixed.panel + mixed.panel:black.def=0")
lieb.main.with.inter.and.lib.hyp

#CREATE TABLE 2 [FIRST TABLE OF REGRESSIONS]
p.table <- apsrtable(
		lieb.main.no.inter, 	lieb.main.with.inter.raw, lieb.main.with.inter, lieb.main.with.inter.and.lib
		#, omitcoef=list ( expression(grep("circuit",coefnames)), expression(grep("year",coefnames)))		
			, omitcoef=list ( expression(grep("circuit",coefnames)))
	, model.names=c("(1)", "(2)", "(3)", "(4)")
	, label="panel_reg_table_1"
	, notes=""
	, caption.position="below"
	  , coef.names=c("Intercept",
	   			 "Panel liberalism",
	   		   "State evidentiary hearing",
	  			 "Out of state representation",
	  			 "Aggravation/Mitigation index",
	  			 "Offender victim index",
	  		 	 "Federal evidentiary hearing",
	  			 "Number of claims",
					 "Year sentenced",
	   		   "Black defendant",
	  			 "Mixed panel",
	  			 "Mixed panel $\\times$ black defendant"
				)
)

p.table #NOTE THAT THE ORDER OF THE COEFFICIENTS DIFFERS FROM THE TABLE IN THE PAPER
#F-tests on B1 + B3=0
c(NA, lieb.main.with.inter.raw.b1b3, lieb.main.with.inter.b1b3, lieb.main.with.inter.and.lib.b1b3)
#p-values:do one-tailed
c(NA, lieb.main.with.inter.raw.hyp[2,"Pr(>Chisq)"]/2  , lieb.main.with.inter.hyp[2,"Pr(>Chisq)"]/2,   lieb.main.with.inter.and.lib.hyp[2,"Pr(>Chisq)"]/2 )

#GETTING PREDICTED PROBABILITIES
nsims <- 1000 #number of simulations

#redo model with Zelig for pred.probs
model.to.use <- zelig(panel.relief ~
                              + mixed.panel	+ black.def	+ mixed.panel:black.def
                            + state.evid.hearing + out.state.representation  + aggrav.mit.index  + offender.victim.index 
                            + fed.evid.hearing  + number.claims +  year.sentence.trend
                            + circuit4 + circuit5 + circuit8 + circuit10 + circuit11 
                            ,	data = lieb.cd,model="logit")
n.cases <- nrow(lieb.cd)
lieb.panel.all.sims <- rep(list( matrix(NA, ncol=n.cases, nrow=nsims )), 4)
names(lieb.panel.all.sims)  <- c("WDWP", "WDBP", "BDWP", "BDBP")
x <- setx(model.to.use,cond=T)
x$black.def <- 0 #white defendant
x$mixed.panel <-0# white panel

set.seed(1223)

#note "WP" = all-white panel= "nonmixed panel' in the paper

lieb.panel.all.sims[["WDWP"]] <- rowMeans( Zelig::sim(model.to.use,x, num=nsims)$qi$ev )

x$black.def <- 0#white defendant,
x$mixed.panel <-1# black panel
lieb.panel.all.sims[["WDBP"]] <- rowMeans(Zelig::sim(model.to.use,x, num=nsims)$qi$ev)
x$black.def <- 1#black defendant,
x$mixed.panel <-0 #white panel
lieb.panel.all.sims[["BDWP"]]  <- rowMeans(Zelig::sim(model.to.use,x, num=nsims)$qi$ev)
x$black.def <- 1#black defendant,
x$mixed.panel <-1# black panel
lieb.panel.all.sims[["BDBP"]]  <-  rowMeans(Zelig::sim(model.to.use,x, num=nsims)$qi$ev)
#get quanitles estimates
lieb.panel.all.estimates <- ldply(lieb.panel.all.sims, quantile, probs=c(.025, .5,.975))
matrix (data= with(lieb.panel.all.estimates, c(lieb.panel.all.estimates[.id=="WDWP",3], lieb.panel.all.estimates[.id=="BDWP",3],
	lieb.panel.all.estimates[.id=="WDBP",3], lieb.panel.all.estimates[.id=="BDBP",3])),	nrow=2, ncol=2, byrow=T,
	dimnames = list(c("White panel", "Mixed panel"),     c("White def.", "Black.def")))
quantile(lieb.panel.all.sims[["BDBP"]]  - lieb.panel.all.sims[["BDWP"]], c(.025,.5,.975) )
mean ( lieb.panel.all.sims[["BDBP"]]  > lieb.panel.all.sims[["BDWP"]] ) #is difference among black defendants signficant?
mean ( lieb.panel.all.sims[["WDBP"]]  > lieb.panel.all.sims[["WDWP"]] )# is difference among white defendants?

#is difference among black defendants larger?
mean ( (lieb.panel.all.sims[["BDBP"]]  - lieb.panel.all.sims[["BDWP"]]) > (lieb.panel.all.sims[["WDBP"]]  - lieb.panel.all.sims[["WDWP"]] ) ) 

#for text, paragraph that starts "Beginning with black defendants"
lieb.panel.all.estimates#

#what is difference for black defendants?
mean(lieb.panel.all.sims[["BDBP"]]  - lieb.panel.all.sims[["BDWP"]])

#what is 95\% estimate of difference for black defendant
quantile  ( lieb.panel.all.sims[["BDBP"]]  - lieb.panel.all.sims[["BDWP"]], c(.025,.975))
#what is 95\% estimate of difference for white defendants
quantile  ( lieb.panel.all.sims[["WDBP"]]  - lieb.panel.all.sims[["WDWP"]], c(.025,.975))


#plot by defendant than panel
pred.means <- lieb.panel.all.estimates[,3]
pred.lower <- lieb.panel.all.estimates[,2]
pred.upper <- lieb.panel.all.estimates[,4]
axis.text <- 1.4
text.size <- 1.5
y.labels <- c( "Non-mixed\npanel","Mixed panel", "Non-mixed\npanel","Mixed panel")
y.offset <- .2
y.offset2 <- .12
point.offset <- .1

#MAKE FIGURE 3
pdf("pred_probs_new.pdf", width=10, height=6)

par(mfrow=c(1,1),mar=c(4,19,1,1))
plot(pred.means, c(1:4)+point.offset, xlim=c(.2,.8),ylim=c(.8,4.2), pch=19,axes=F, main="", xlab="", ylab="", cex=1.5)
segments(pred.lower, c(1:4)+point.offset,pred.upper, c(1:4)+point.offset, lwd=2)

axis(1,at=seq(.2,.8,.1), labels=c(".2",".3",".4",".5",".6",".7",".8"), mgp=c(2,.5,0), cex.axis=axis.text)
axis(2, at=c(1:4), cex.axis=axis.text, labels=y.labels,las=1, cex.axis=axis.text)
text(par()$usr[1]-y.offset, 3.5, '{', srt = 0, cex = 12,lwd=.5, xpd=NA)
text(par()$usr[1]-y.offset-y.offset2,3.5,"Black\ndefendants",cex=text.size,xpd=NA,font=3)
text(par()$usr[1]-y.offset, y = 1.5, '{', srt = 0, cex = 15 ,xpd=NA)
text(par()$usr[1]-y.offset-y.offset2,1.5,"White\ndefendants",cex=text.size,xpd=NA,font=3)
mtext("Probability of panel relief", 1, line=2.5, cex=2)

dev.off()


#RULING OUT ALTERNATIVE IDEOLOGICAL STORY
#using panel liberalism and panel indicators
use <- subset(lieb.cd, mixed.panel==0)	#DROP MIXED PANELS TO ISOLATE EFFECT OF IDEOLOGY
lieb.robust.liberals.1 <- glm(panel.relief ~
			+ state.evid.hearing + out.state.representation  + aggrav.mit.index  + offender.victim.index 
				+ fed.evid.hearing  + number.claims + year.sentence.trend
				+ circuit4 + circuit5 + circuit8 + circuit10 + circuit11
				+ panel.liberalism*black.def ,data=use,family=binomial(link="logit"))
summary(lieb.robust.liberals.1)

lieb.robust.liberals.2 <- glm(panel.relief ~
			+ state.evid.hearing + out.state.representation  + aggrav.mit.index  + offender.victim.index 
				+ fed.evid.hearing  + number.claims + year.sentence.trend
				+ circuit4 + circuit5 + circuit8 + circuit10 + circuit11
				+ black.def +  rrd.panel + ddr.panel + ddd.panel
				+ rrd.panel:black.def  + ddr.panel:black.def  + ddd.panel:black.def ,data=use,family=binomial(link="logit"))
summary(lieb.robust.liberals.2)

#NOTE THAT THE ORDER OF THE COEFFICIENTS DIFFERS FROM THE TABLE IN THE PAPER
apsrtable(lieb.robust.liberals.1, lieb.robust.liberals.2
	, model.names=c("(1)","(2)")
	, label="ideological_reg_table"
	, notes=""
	, caption.position="below"
	, omitcoef=list ( expression(grep("circuit",coefnames)))		
	, coef.names=c("Intercept", 
		"State evidentiary hearing", 
		"Out of state representation", 	
		"Aggravation/Mitigation index", 
		 "Offender victim index", 
		 "Federal evidentiary hearing", 
		 "Number of claims", 
			"Year sentenced",
			"Panel liberalism",
			"Black defendant", 
			"Panel liberalism $\\times$ black defendant", 
			"RRD panel",
			"DDR panel", 	 
			"DDD panel",
    	"RRD panel $\\times$ black defendant",
			"DDR panel  $\\times$ black defendant",
			"DDD panel  $\\times$ black defendant"),
			caption="Regression models"
		)

######################################################################################
######################################################################################
#ROBUSTNESS CHECKS in Appendix
######################################################################################
#Now cross-validate dropping judges 1 by 1
set.seed(153)

temp.data <- NA
temp.coefs <- NA

use <- subset(lieb.cd)
x<-with(use, table(code1[black1==1]))
y<-with(use, table(code2[black2==1]))
z<-with(use, table(code3[black3==1]))

judges.to.exclude.lieb <- sort(as.numeric(unique(c(names(x),names(y),names(z)))))
judges.names <- c("Higginbotham", "Hatchett", "McMillian", "Farris", "Poole", "Davis")

#loop over each judge
for (i in 1:6){
	temp.data <- subset(lieb.cd, code1!=judges.to.exclude.lieb[i] & code2!=judges.to.exclude.lieb[i] & code3!=judges.to.exclude.lieb[i])
	tempname <- judges.names[i]
 	model.to.use <-	assign(paste("lieb.cv.panel.model.no.party.",i,sep=""), zelig(panel.relief ~
			+ black.def
			+ mixed.panel
			+ mixed.panel:black.def  
			+ state.evid.hearing + out.state.representation  + aggrav.mit.index  + offender.victim.index 
			+ fed.evid.hearing  + number.claims + year.sentence.trend
			+ circuit4 + circuit5 + circuit8 + circuit10 + circuit11
	 		, data = temp.data, model="logit") )
#f-test
assign(paste("lieb.f.test.black.def.",i,sep=""), linearHypothesis(model.to.use, "mixed.panel+ black.def:mixed.panel =0"))
	#do sims in loop/get pred probs
	assign(paste("lieb.white.def.white.panel.sims.",i,sep=""), matrix(NA, ncol=length(model.to.use$y), nrow=nsims ))
	assign(paste("lieb.white.def.black.panel.sims.",i,sep=""), matrix(NA, ncol=length(model.to.use$y), nrow=nsims ))	
	assign(paste("lieb.black.def.white.panel.sims.",i,sep=""), matrix(NA, ncol=length(model.to.use$y), nrow=nsims ))
	assign(paste("lieb.black.def.black.panel.sims.",i,sep=""), matrix(NA, ncol=length(model.to.use$y), nrow=nsims ))
	x <- setx(model.to.use,cond=T)
	#white defendant, white panel
	temp.x <- x
	temp.x$black.def <- 0
	temp.x$mixed.panel <-0
	sim.temp <- Zelig::sim(model.to.use,temp.x, num=nsims)
	preds <- sim.temp$qi$ev
	assign(paste("lieb.white.def.white.panel.sims.",i,sep=""),rowMeans(preds))
	temp <- 	assign(paste("lieb.white.def.white.panel.sims.",i,sep=""),rowMeans(preds))
	assign(paste("lieb.white.def.white.panel.estimates.",i,sep=""), quantile(temp, probs=c(.025, .5,.975)))
	#white defendant, black panel
	temp.x <- x
	temp.x$black.def <- 0
	temp.x$mixed.panel <-1
	sim.temp <- Zelig::sim(model.to.use,temp.x, num=nsims)
	preds <- sim.temp$qi$ev
	assign(paste("lieb.white.def.black.panel.sims.",i,sep=""),rowMeans(preds))
	temp <- 	assign(paste("lieb.white.def.black.panel.sims.",i,sep=""),rowMeans(preds))
	assign(paste("lieb.white.def.black.panel.estimates.",i,sep=""), quantile(temp, probs=c(.025, .5,.975)))
	#black defendant, white panel
	temp.x <- x
	temp.x$black.def <- 1
	temp.x$mixed.panel <-0
	sim.temp <- Zelig::sim(model.to.use,temp.x, num=nsims)
	preds <- sim.temp$qi$ev
	assign(paste("lieb.black.def.white.panel.sims.",i,sep=""),rowMeans(preds))
	temp <- 	assign(paste("lieb.black.def.white.panel.sims.",i,sep=""),rowMeans(preds))
	assign(paste("lieb.black.def.white.panel.estimates.",i,sep=""), quantile(temp, probs=c(.025, .5,.975)))	
	#black defendant, black panel
	temp.x <- x
	temp.x$black.def <- 1
	temp.x$mixed.panel <-1
	sim.temp <- Zelig::sim(model.to.use,temp.x, num=nsims)
	preds <- sim.temp$qi$ev
	assign(paste("lieb.black.def.black.panel.sims.",i,sep=""),rowMeans(preds))
	temp <- 	assign(paste("lieb.black.def.black.panel.sims.",i,sep=""),rowMeans(preds))
	assign(paste("lieb.black.def.black.panel.estimates.",i,sep=""), quantile(temp, probs=c(.025, .5,.975)))	
		}
		
#make regression table	#NOTE THAT THE ORDER OF THE COEFFICIENTS DIFFERS FROM THE TABLE IN THE PAPER	
#################
panel.cv.table.lieb <- apsrtable(		lieb.cv.panel.model.no.party.1, lieb.cv.panel.model.no.party.2,
		lieb.cv.panel.model.no.party.3,lieb.cv.panel.model.no.party.4,lieb.cv.panel.model.no.party.5,lieb.cv.panel.model.no.party.6
		, model.names=c("(1)","(2)","(3)","(4)", "(5)", "(6)")
		, label="panel_reg_table_cv_lieb"
		, notes=""
		, caption.position="below"
		, omitcoef= c("circuit4", "circuit5", "circuit8", "circuit10", "circuit11", "circuit9")
		, coef.names=c("Intercept", "Black defendant", "Mixed panel", 
				"State evidentiary hearing", 
				"Out of state representation", 	"Aggravation/Mitigation index", 
				"Offender victim index", "Federal evidentiary hearing", "Number of claims", "Year sentenced", 
				"Mixed panel $\\times$ black defendant")
		, caption="Models of panel decisions to grant relief using the Liebman data.	In each model, one of the six unique black judges in the data is excluded.
		Fixed effects	for circuits estimated in each model but not displayed.
	 	$^*$ indicates significance at $p< 0.05$."
		)
panel.cv.table.lieb
#add f-tests by hand
c(lieb.f.test.black.def.1$"Pr(>Chisq)"[2], lieb.f.test.black.def.2$"Pr(>Chisq)"[2],
	lieb.f.test.black.def.3$"Pr(>Chisq)"[2],lieb.f.test.black.def.4$"Pr(>Chisq)"[2],
	lieb.f.test.black.def.5$"Pr(>Chisq)"[2],lieb.f.test.black.def.6$"Pr(>Chisq)"[2])/2#one-tailed
	
#plot each prediction

pdf("pred_probs_cv.pdf", width=8.5, height=8)

par(mfrow=c(1,1),mar=c(4,19,1,1))
y <- c(1:4)
offset <- seq(-.4,.4,length.out=13)
point.type <-19
plot(y, y, ylim=c(.5,4.5), xlim=c(.2,.8), type="n", axes=F,  main="", xlab="", ylab="", cex=1.5)

#lieb
for (i in 1:6){
	#White def white panel
	use <- get(paste0("lieb.white.def.white.panel.estimates.",i))
	points(use[2], y[1]+offset[i], pch=point.type)
	segments(use[1], y[1]+offset[i],use[3], y[1]+offset[i], lwd=2)
	#White def black panel
	use <- get(paste0("lieb.white.def.black.panel.estimates.",i))
	points(use[2], y[2]+offset[i], pch=point.type)
	segments(use[1], y[2]+offset[i],use[3], y[2]+offset[i], lwd=2)
	#Black def white panel
	use <- get(paste0("lieb.black.def.white.panel.estimates.",i))
	points(use[2], y[3]+offset[i], pch=point.type)
	segments(use[1], y[3]+offset[i],use[3], y[3]+offset[i], lwd=2)
	#Black def black panel
	use <- get(paste0("lieb.black.def.black.panel.estimates.",i))
	points(use[2], y[4]+offset[i], pch=point.type)
	segments(use[1], y[4]+offset[i],use[3], y[4]+offset[i], lwd=2)	
	}

y.offset <- .28
y.offset2 <- .2
axis(1,at=seq(.2,.8,.1), labels=c(".2",".3",".4",".5",".6",".7",".8"), mgp=c(2,.5,0), cex.axis=axis.text)
axis(2, at=c(1:4), cex.axis=axis.text, labels=y.labels,las=1, cex.axis=axis.text)
text(par()$usr[1]-y.offset, 3.5, '{', srt = 0, cex = 12,lwd=.5, xpd=NA)
text(par()$usr[1]-y.offset-y.offset2+.03,3.5,"Black\ndefendants",cex=text.size,xpd=NA,font=3)
text(par()$usr[1]-y.offset, y = 1.5, '{', srt = 0, cex = 15 ,xpd=NA)
text(par()$usr[1]-y.offset-y.offset2+.03,1.5,"White\ndefendants",cex=text.size,xpd=NA,font=3)
mtext("Probability of panel relief", 1, line=2.5, cex=2)

dev.off()

# ####	
# In addition, even when using Model (3)---which contains the noisiest estimates of $\beta_1+\beta_3$---in 90\% of simulations,
# the probability of relief when the panel is mixed (conditional on the defendant being black) is higher than the probability of 
# relief when the panel is non-mixed.

mean(lieb.black.def.black.panel.sims.3 > lieb.black.def.white.panel.sims.3)

